## test out actual values with differing sigR, # of stocks, sample size
dep <- .2

mu <- 0
tau <- 1

T <- 40
m <- 10

L <- 4
qI <- 1; sigI <- 0.5; sigN <- 0.5; qN <- 1

pars <- list(mu=mu, tau=tau, qI=qI, sigI=sigI, sigN=sigN, qN=qN)

sims <- 10

system("admb -r simulate_RE"); system("admb -r estimate_nlme");

meso.results <- list(); I <- 1
for (m in c(10, 20)) {
  pars$m <- m

  B0 <- rep(40, m); p <- rep(1, m); w <- rep(1, m); M <- rep(0.2, m); q <- rep(1, m)

  phi0 <- rep(0,m) # mle estimates
  for (k in 1:m)
    phi0[k] <- (1 - p[k]*w[k]*exp(-M[k])) /
      (1 - (1+p[k])*exp(-M[k]) + p[k] * exp(-2*M[k]))

  pars$B0 <- B0; pars$p <- p; pars$w <- w; pars$M <- M; pars$q <- q; pars$phi0 <- phi0
  
  for (T in c(40)) {
    pars$T <- T
    T_dep <- T/2

    Eb <- with(pars, rep(c(rep(0, len=T/8), seq(0,.1, len=T/8), seq(.1,.25, len=T/4), seq(.25, 0, len=T/4), rep(0, len=T/4), 0)))
    E <- with(pars, matrix(rep(Eb, times = m), T+1, m, byrow = F))
    pars$E <- E

    results <- list()
    for (sigR in c(.2)) {
      pars$sigR <- sigR
      
      for (i in 1:sims) {
        pars$beta <- rnorm(m, mu, tau)
        sim <- with(pars, run_opmodel(pars, m, T, E, L, phi0))
        out <- run_meta(est, sim, pars, method = "nlme", real = TRUE)

        results[[i]] <- data.frame(beta = pars$beta,
                                   beta_est = with(out, eps * exp(log_tau) + mu),
                                   eps = out$eps,
                                   B0 = B0,
                                   B0_est = out$B0,
                                   mu = mu,
                                   mu_est = out$mu,
                                   tau = tau,
                                   tau_est= exp(out$log_tau),
                                   sigR = sigR,
                                   ##sigR_sim = as.numeric(lapply(est, function(x) x$sigR)),
                                   sigR_est = exp(out$log_sigR))
      }

      meso.results[[I]] <- results
      I <- I + 1
    }
  }
}

## plot
aho <- do.call(rbind, meso.results[[1]])
baka <- do.call(rbind, meso.results[[2]])

par(mfrow = c(2,1))
plot(aho$beta, aho$beta_est)
plot(baka$beta, baka$beta_est)

plot(unique(aho$mu_est), unique(aho$tau_est))
plot(unique(baka$mu_est), unique(baka$tau_est))

dotchart(unique(aho$mu_est), xlim = c(-.5,.5)); abline(v = aho$mu)
dotchart(unique(baka$mu_est), xlim = c(-.5,.5)); abline(v = baka$mu)

dotchart(unique(aho$tau_est), xlim = 1+c(-.5,.5)); abline(v = aho$tau)
dotchart(unique(baka$tau_est), xlim = 1+c(-.5,.5)); abline(v = baka$tau)

dotchart(aho$beta_est - aho$beta); abline(v = 0)
dotchart(baka$beta_est - baka$beta); abline(v = 0)

dotchart(aho$sigR_est - aho$sigR, xlim = c(-.07,.07)); abline(v = 0)
dotchart(baka$sigR_est - baka$sigR, xlim = c(-.07,.07)); abline(v = 0)

boxplot(aho$beta_est - aho$beta ~ as.character(aho$mu_est))
boxplot(baka$beta_est - baka$beta ~ as.character(baka$mu_est))
